function stop = h2out(tt, uu, flag, params)
  stop = 0;
  
  if (~params.plotHydro)
    return
  end

  global h2out_lasttime h2out_points h2out_lastu h2out_nlines;
  if strcmp(flag, 'init')
    h2out_lastu = zeros(2, params.Nx, params.Ny);
    h2out_lasttime = 0;
    h2out_nlines = 12;
    h2out_points = zeros(h2out_nlines*(params.Nx+params.Ny), 2);

    for l=1:h2out_nlines
      h2out_points(h2out_nlines*params.Nx + (l-1)*params.Ny+1 : h2out_nlines*params.Nx + l*params.Ny, 1) = (l-1)*params.Lx/h2out_nlines;
      h2out_points(h2out_nlines*params.Nx + (l-1)*params.Ny+1 : h2out_nlines*params.Nx + l*params.Ny, 2) = (0:params.Ny-1)/params.Ny*params.Ly;

      h2out_points((l-1)*params.Nx+1 : l*params.Nx, 1) = (0:params.Nx-1)/params.Nx*params.Lx;
      h2out_points((l-1)*params.Nx+1 : l*params.Nx, 2) = (l-1)*params.Ly/h2out_nlines;
    end
  end


  
  if strcmp(flag, 'done')
    return
  end

  Nx = params.Nx;
  Ny = params.Ny;
  Nt = Nx.*Ny;

  tsteps = size(uu,2);
  disp(sprintf('t = %f, Nt = %d', tt(tsteps), tsteps));

  u = uu(:, tsteps);
  t = tt(:, tsteps);


 
  T  = reshape(u(0*Nt+1 : 1*Nt), Nx, Ny);
  ux = reshape(u(1*Nt+1 : 2*Nt), Nx, Ny);
  uy = reshape(u(2*Nt+1 : 3*Nt), Nx, Ny);
  P1 = reshape(u(3*Nt+1 : 4*Nt), Nx, Ny);
  P2 = reshape(u(4*Nt+1 : 5*Nt), Nx, Ny);

  Tdx  = dx(T);
  Tdy  = dy(T);
  uxdx = dx(ux);
  uxdy = dy(ux);
  uydx = dx(uy);
  uydy = dy(uy);
  uxdxdx  = dx(uxdx);
  uxdxdy  = dy(uxdx);
  uxdydy  = dy(uxdy);

  uydxdx  = dx(uydx);
  uydxdy  = dy(uydx);
  uydydy  = dy(uydy);

  T1 = sqrt( (Tdx.^2+Tdy.^2) )./T;

  u2 = 0;
  u1 = sqrt(uxdx.^2 + uxdy.^2);
  u2 = sqrt(uydx.^2 + uydy.^2);
  %u1 = ((ux.^2+uy.^2).*(1+ux.^2+uy.^2).^(-1)).^(-1/2).*((ux.^2+uy.^2).^(-1).*(1+ux.^2+uy.^2).^(-3).*(ux.^2.*(uxdx.^2+uxdy.^2)+2.*ux.*uy.*(uxdx.*uydx+uxdy.*uydy)+uy.^2.*(uydx.^2+uydy.^2))).^(1/2);
  %u2 = (1/2).*((ux.^2+uy.^2).*(1+ux.^2+uy.^2).^(-1)).^(1/2).*((ux.^2+uy.^2).^(-1).*(1+ux.^2+uy.^2).^(-3).*(ux.^2.*(uxdx.^2+uxdy.^2)+2.*ux.*uy.*(uxdx.*uydx+uxdy.*uydy)+uy.^2.*(uydx.^2+uydy.^2))).^(-1/2).*((ux.^2+uy.^2).^(-4).*(1+ux.^2+uy.^2).^(-4).*(ux.^2.*(uxdx.^2+uxdy.^2)+2.*ux.*uy.*(uxdx.*uydx+uxdy.*uydy)+uy.^2.*(uydx.^2+uydy.^2)).^(-1).*((2.*(ux.*uxdx+uy.*uydx).*(ux.^2.*(uxdx.^2+uxdy.^2)+2.*ux.*uy.*(uxdx.*uydx+uxdy.*uydy)+uy.^2.*(uydx.^2+uydy.^2))+6.*(ux.^2+uy.^2).*(ux.*uxdx+uy.*uydx).*(ux.^2.*(uxdx.^2+uxdy.^2)+2.*ux.*uy.*(uxdx.*uydx+uxdy.*uydy)+uy.^2.*(uydx.^2+uydy.^2))+2.*(1+ux.^2+uy.^2).*(ux.*uxdx+uy.*uydx).*(ux.^2.*(uxdx.^2+uxdy.^2)+2.*ux.*uy.*(uxdx.*uydx+uxdy.*uydy)+uy.^2.*(uydx.^2+uydy.^2))+(-2).*(ux.^2+uy.^2).*(1+ux.^2+uy.^2).*(ux.^2.*(uxdx.*uxdxdx+uxdxdy.*uxdy)+uy.*(uydx.*(uxdx.^2+uydx.^2+uy.*uydxdx)+uxdx.*uxdy.*uydy+uy.*uydxdy.*uydy+uydx.*uydy.^2)+ux.*(uxdx.^3+uxdx.*uxdy.^2+uxdx.*(uydx.^2+uy.*uydxdx)+uy.*(uxdxdx.*uydx+uxdxdy.*uydy)+uxdy.*(uy.*uydxdy+uydx.*uydy)))).^2+(2.*(ux.*uxdy+uy.*uydy).*(ux.^2.*(uxdx.^2+uxdy.^2)+2.*ux.*uy.*(uxdx.*uydx+uxdy.*uydy)+uy.^2.*(uydx.^2+uydy.^2))+6.*(ux.^2+uy.^2).*(ux.*uxdy+uy.*uydy).*(ux.^2.*(uxdx.^2+uxdy.^2)+2.*ux.*uy.*(uxdx.*uydx+uxdy.*uydy)+uy.^2.*(uydx.^2+uydy.^2))+2.*(1+ux.^2+uy.^2).*(ux.*uxdy+uy.*uydy).*(ux.^2.*(uxdx.^2+uxdy.^2)+2.*ux.*uy.*(uxdx.*uydx+uxdy.*uydy)+uy.^2.*(uydx.^2+uydy.^2))+(-2).*(ux.^2+uy.^2).*(1+ux.^2+uy.^2).*(ux.^2.*(uxdx.*uxdxdy+uxdy.*uxdydy)+uy.*(uxdx.*uxdy.*uydx+uy.*uydx.*uydxdy+uxdy.^2.*uydy+uydy.^3+uydy.*(uydx.^2+uy.*uydydy))+ux.*(uxdy.^3+uxdx.*uydx.*uydy+uy.*(uxdxdy.*uydx+uxdx.*uydxdy+uxdydy.*uydy)+uxdy.*(uxdx.^2+uydy.^2+uy.*uydydy)))).^2)).^(1/2);





  disp(sprintf('hydro-check: T: (%.0e, %.0e), u: (%.0e, %.0e), u2: (%.0e, %.0e)', mean(mean(T1)), max(max(T1)), mean(mean(u1)), max(max(u1)), mean(mean(u2)), max(max(u2)) ));

  w = uxdy - uydx;
  w = w * 5;

  if (false)
    dt = t - h2out_lasttime;
    h2out_lasttime = t;
    for n=1:h2out_nlines*(params.Nx+params.Ny)
      lx = round(h2out_points(n,1)/params.Lx*params.Nx) + 1;
      ly = round(h2out_points(n,2)/params.Ly*params.Ny) + 1;
      if lx > params.Nx
        lx = lx - params.Nx;
      end
      if ly > params.Ny
        ly = ly - params.Ny;
      end
      if lx < 1
        lx = lx + params.Nx;
      end
      if ly < 1
        ly = ly + params.Ny;
      end
      h2out_points(n,1) = h2out_points(n,1) + h2out_lastu(1,lx,ly) .* dt;
      h2out_points(n,2) = h2out_points(n,2) + h2out_lastu(2,lx,ly) .* dt;
    end
    ut = sqrt(1 + ux.^2 + uy.^2);
    h2out_lastu(1,:,:) = ux./ut;
    h2out_lastu(2,:,:) = uy./ut;
  end




  figure(1);

  plotsx = 3;
  plotsy = 2;

  f = 1;

  subplot(plotsy, plotsx, 1);
  T = T - params.TempOffset;
  maxT = max(max(abs(T)));
  meanT = mean(mean(T));
  T = T ./ maxT;
  image(T.*params.colormapOffset*f+params.colormapOffset); 
  axis square;
  title(sprintf('T  @  t = %.1f, DT = %.5f', t, meanT));
  if false
  hold on
  gridskip = 4;
  [x, y] = meshgrid(1:gridskip:params.Nx,1:gridskip:params.Ny);
  gridux = ux(1:gridskip:params.Nx, 1:gridskip:params.Ny);
  griduy = uy(1:gridskip:params.Nx, 1:gridskip:params.Ny);
  quiver(x, y, griduy, gridux);
  axis([1 params.Nx 1 params.Ny]);
  hold off
  end




  fv = 1;
  subplot(plotsy, plotsx, 2);
  maxux = max(max(abs(ux)));
  image(ux./maxux*params.colormapOffset*fv +params.colormapOffset);
  axis square;
  title(sprintf('ux,  max  = %.4f', maxux));
  % if (VIDEO) & (frame_active)
  %   frame = getframe;
  %   writeVideo(vvx,frame);
  % end

  subplot(plotsy, plotsx, 3);
  maxuy = max(max(abs(uy)));
  image(uy./maxuy*params.colormapOffset*fv +params.colormapOffset);
  axis square;
  title(sprintf('uy,  max  = %.4f', maxuy));
  % if (VIDEO) & (frame_active)
  %   frame = getframe;
  %   writeVideo(vvy,frame);
  % end

  subplot(plotsy, plotsx, 4);
  image(w*f*params.colormapOffset +params.colormapOffset);
  axis square;
  title('w');
  % if (VIDEO) & (frame_active)
  %   frame = getframe;
  %   writeVideo(vvy,frame);
  % end

  subplot(plotsy, plotsx, 5);
  image(P1*params.colormapOffset/8+params.colormapOffset);
  axis square;
  title('P1');
  if (false)
    plot(h2out_points(:,1), h2out_points(:,2), '.', 'MarkerSize', 1);
    axis([0 params.Lx 0 params.Ly]);
  end
  maxw = max(max(abs(w)));
  w = w/maxw*params.colormapOffset + params.colormapOffset;
  image(w);
  axis square;
  title('w auto');

  maxp2 = max(max(abs(P2)));
  subplot(plotsy, plotsx, 6);
  image(P2/maxp2*params.colormapOffset+params.colormapOffset);
  axis square;
  title('P2');


  refresh();
end
